home *** CD-ROM | disk | FTP | other *** search
/ Aminet 3 / Aminet 3 - July 1994.iso / Aminet / mus / misc / mpegaudio.lha / mpegaudio / psy.c < prev    next >
C/C++ Source or Header  |  1994-03-21  |  19KB  |  444 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. psy.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@3d.enet.dec.com             *
  10.  *                                                                    *
  11.  * VERSION 3.9                                                        *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Davis Pan           start of version 1.0 records          *
  15.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  16.  * 7/10/91  Earle Jennings      Ported to MsDos.                      *
  17.  *                              replace of floats with FLOAT          *
  18.  * 2/11/92  W. Joseph Carter    Fixed mem_alloc() arg for "absthr".   *
  19.  * 7/24/92  M. Iwadare          HANN window coefficients modified.    *
  20.  * 7/27/92  Masahiro Iwadare    Bug fix, FFT modification for Layer 3 *
  21.  * 7/27/92  Masahiro Iwadare    Bug fix, "new", "old", and "oldest"   *
  22.  *                              updates                               *
  23.  * 8/07/92  Mike Coleman        Bug fix, read_absthr()                *
  24.  **********************************************************************/
  25.  
  26. #include "common.h"
  27. #include "encoder.h"
  28.  
  29. void psycho_anal(buffer,savebuf,chn,lay,snr32,sfreq)
  30. short int *buffer;
  31. short int savebuf[1056];
  32. int   chn, lay;
  33. FLOAT snr32[32];
  34. double sfreq;        /* to match prototype : float args are always double */
  35. {
  36.  unsigned int   i, j, k;
  37.  FLOAT          r_prime, phi_prime;
  38.  FLOAT          freq_mult, bval_lo, minthres, sum_energy;
  39.  double         tb, temp1, temp2, temp3;
  40.  
  41. /* The static variables "r", "phi_sav", "new", "old" and "oldest" have    */
  42. /* to be remembered for the unpredictability measure.  For "r" and        */
  43. /* "phi_sav", the first index from the left is the channel select and     */
  44. /* the second index is the "age" of the data.                             */
  45.  
  46.  static int     new = 0, old = 1, oldest = 0;
  47.  static int     init = 0, flush, sync_flush, syncsize, sfreq_idx;
  48.  
  49. /* The following static variables are constants.                           */
  50.  
  51.  static double  nmt = 5.5;
  52.  
  53.  static FLOAT   crit_band[27] = {0,  100,  200, 300, 400, 510, 630,  770,
  54.                                920, 1080, 1270,1480,1720,2000,2320, 2700,
  55.                               3150, 3700, 4400,5300,6400,7700,9500,12000,
  56.                              15500,25000,30000};
  57.  
  58.  static FLOAT   bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
  59.                             10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,
  60.                              4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,
  61.                              4.5,  4.5,  4.5,  3.5,  3.5,  3.5};
  62.  
  63. /* The following pointer variables point to large areas of memory         */
  64. /* dynamically allocated by the mem_alloc() function.  Dynamic memory     */
  65. /* allocation is used in order to avoid stack frame or data area          */
  66. /* overflow errors that otherwise would have occurred at compile time     */
  67. /* on the Macintosh computer.                                             */
  68.  
  69.  FLOAT          *grouped_c, *grouped_e, *nb, *cb, *ecb, *bc;
  70.  FLOAT          *wsamp_r, *wsamp_i, *phi, *energy;
  71.  FLOAT          *c, *fthr;
  72.  F32            *snrtmp;
  73.  
  74.  static int     *numlines;
  75.  static int     *partition;
  76.  static FLOAT   *cbval, *rnorm;
  77.  static FLOAT   *window;
  78.  static FLOAT   *absthr;
  79.  static double  *tmn;
  80.  static FCB     *s;
  81.  static FHBLK   *lthr;
  82.  static F2HBLK  *r, *phi_sav;
  83.  
  84. /* These dynamic memory allocations simulate "automatic" variables        */
  85. /* placed on the stack.  For each mem_alloc() call here, there must be    */
  86. /* a corresponding mem_free() call at the end of this function.           */
  87.  
  88.  grouped_c = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_c");
  89.  grouped_e = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_e");
  90.  nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb");
  91.  cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb");
  92.  ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb");
  93.  bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc");
  94.  wsamp_r = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_r");
  95.  wsamp_i = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_i");
  96.  phi = (FLOAT *) mem_alloc(sizeof(FBLK), "phi");
  97.  energy = (FLOAT *) mem_alloc(sizeof(FBLK), "energy");
  98.  c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c");
  99.  fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr");
  100.  snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");
  101.  
  102.  if(init==0){
  103.  
  104. /* These dynamic memory allocations simulate "static" variables placed    */
  105. /* in the data space.  Each mem_alloc() call here occurs only once at     */
  106. /* initialization time.  The mem_free() function must not be called.      */
  107.  
  108.      numlines = (int *) mem_alloc(sizeof(ICB), "numlines");
  109.      partition = (int *) mem_alloc(sizeof(IHBLK), "partition");
  110.      cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
  111.      rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
  112.      window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");
  113.      absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr");
  114.      tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
  115.      s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
  116.      lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
  117.      r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");
  118.      phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");
  119.  
  120.      i = sfreq + 0.5;
  121.      switch(i){
  122.         case 32000: sfreq_idx = 0; break;
  123.         case 44100: sfreq_idx = 1; break;
  124.         case 48000: sfreq_idx = 2; break;
  125.         default:    printf("error, invalid sampling frequency: %d Hz\n",i);
  126.         exit(-1);
  127.      }
  128.      printf("absthr[][] sampling frequency index: %d\n",sfreq_idx);
  129.      read_absthr(absthr, sfreq_idx);
  130.      if(lay==1){
  131.         flush = 384;
  132.         syncsize = 1024;
  133.         sync_flush = 576;
  134.      }
  135.      else {
  136.         flush = 384*3.0/2.0;
  137.         syncsize = 1056;
  138.         sync_flush = syncsize - flush;
  139.      }
  140. /* calculate HANN window coefficients */
  141. /*   for(i=0;i<BLKSIZE;i++)window[i]=0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0))); */
  142.      for(i=0;i<BLKSIZE;i++)window[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));
  143. /* reset states used in unpredictability measure */
  144.      for(i=0;i<HBLKSIZE;i++){
  145.         r[0][0][i]=r[1][0][i]=r[0][1][i]=r[1][1][i]=0;
  146.         phi_sav[0][0][i]=phi_sav[1][0][i]=0;
  147.         phi_sav[0][1][i]=phi_sav[1][1][i]=0;
  148.         lthr[0][i] = 60802371420160.0;
  149.         lthr[1][i] = 60802371420160.0;
  150.      }
  151. /*****************************************************************************
  152.  * Initialization: Compute the following constants for use later             *
  153.  *    partition[HBLKSIZE] = the partition number associated with each        *
  154.  *                          frequency line                                   *
  155.  *    cbval[CBANDS]       = the center (average) bark value of each          *
  156.  *                          partition                                        *
  157.  *    numlines[CBANDS]    = the number of frequency lines in each partition  *
  158.  *    tmn[CBANDS]         = tone masking noise                               *
  159.  *****************************************************************************/
  160. /* compute fft frequency multiplicand */
  161.      freq_mult = sfreq/BLKSIZE;
  162.  
  163. /* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/
  164.      for(i=0;i<HBLKSIZE;i++){
  165.         temp1 = i*freq_mult;
  166.         j = 1;
  167.         while(temp1>crit_band[j])j++;
  168.         fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);
  169.      }
  170.      partition[0] = 0;
  171. /* temp2 is the counter of the number of frequency lines in each partition */
  172.      temp2 = 1;
  173.      cbval[0]=fthr[0];
  174.      bval_lo=fthr[0];
  175.      for(i=1;i<HBLKSIZE;i++){
  176.         if((fthr[i]-bval_lo)>0.33){
  177.            partition[i]=partition[i-1]+1;
  178.            cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  179.            cbval[partition[i]] = fthr[i];
  180.            bval_lo = fthr[i];
  181.            numlines[partition[i-1]] = temp2;
  182.            temp2 = 1;
  183.         }
  184.         else {
  185.            partition[i]=partition[i-1];
  186.            cbval[partition[i]] += fthr[i];
  187.            temp2++;
  188.         }
  189.      }
  190.      numlines[partition[i-1]] = temp2;
  191.      cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
  192.  
  193. /************************************************************************
  194.  * Now compute the spreading function, s[j][i], the value of the spread-*
  195.  * ing function, centered at band j, for band i, store for later use    *
  196.  ************************************************************************/
  197.      for(j=0;j<CBANDS;j++){
  198.         for(i=0;i<CBANDS;i++){
  199.            temp1 = (cbval[i] - cbval[j])*1.05;
  200.            if(temp1>=0.5 && temp1<=2.5){
  201.               temp2 = temp1 - 0.5;
  202.               temp2 = 8.0 * (temp2*temp2 - 2.0 * temp2);
  203.            }
  204.            else temp2 = 0;
  205.            temp1 += 0.474;
  206.            temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));
  207.            if(temp3 <= -100) s[i][j] = 0;
  208.            else {
  209.               temp3 = (temp2 + temp3)*LN_TO_LOG10;
  210.               s[i][j] = exp(temp3);
  211.            }
  212.         }
  213.      }
  214.  
  215.   /* Calculate Tone Masking Noise values */
  216.      for(j=0;j<CBANDS;j++){
  217.         temp1 = 15.5 + cbval[j];
  218.         tmn[j] = (temp1>24.5) ? temp1 : 24.5;
  219.   /* Calculate normalization factors for the net spreading functions */
  220.         rnorm[j] = 0;
  221.         for(i=0;i<CBANDS;i++){
  222.            rnorm[j] += s[j][i];
  223.         }
  224.      }
  225.      init++;
  226.  }
  227.  
  228. /************************* End of Initialization *****************************/
  229.  switch(lay) {
  230.   case 1:
  231.   case 2:
  232.      for(i=0; i<lay; i++){
  233. /*****************************************************************************
  234.  * Net offset is 480 samples (1056-576) for layer 2; this is because one must*
  235.  * stagger input data by 256 samples to synchronize psychoacoustic model with*
  236.  * filter bank outputs, then stagger so that center of 1024 FFT window lines *
  237.  * up with center of 576 "new" audio samples.                                *
  238.  *                                                                           *
  239.  * For layer 1, the input data still needs to be staggered by 256 samples,   *
  240.  * then it must be staggered again so that the 384 "new" samples are centered*
  241.  * in the 1024 FFT window.  The net offset is then 576 and you need 448 "new"*
  242.  * samples for each iteration to keep the 384 samples of interest centered   *
  243.  *****************************************************************************/
  244.         for(j=0; j<syncsize; j++){
  245.            if(j<(sync_flush))savebuf[j] = savebuf[j+flush];
  246.            else savebuf[j] = *buffer++;
  247.            if(j<BLKSIZE){
  248. /**window data with HANN window***********************************************/
  249.               wsamp_r[j] = window[j]*((FLOAT) savebuf[j]);
  250.               wsamp_i[j] = 0;
  251.            }
  252.         }
  253. /**Compute FFT****************************************************************/
  254.         fft(wsamp_r,wsamp_i,energy,phi,1024);
  255. /*****************************************************************************
  256.  * calculate the unpredictability measure, given energy[f] and phi[f]        *
  257.  *****************************************************************************/
  258. /*only update data "age" pointers after you are done with both channels      */
  259. /*for layer 1 computations, for the layer 2 double computations, the pointers*/
  260. /*are reset automatically on the second pass                                 */
  261.          if(lay==2 || (lay==1 && chn==0) ){
  262.            if(new==0){new = 1; oldest = 1;}
  263.            else {new = 0; oldest = 0;}
  264.            if(old==0)old = 1; else old = 0;
  265.         }
  266.         for(j=0; j<HBLKSIZE; j++){
  267.            r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];
  268.            phi_prime = 2.0 * phi_sav[chn][old][j] - phi_sav[chn][oldest][j];
  269.            r[chn][new][j] = sqrt((double) energy[j]);
  270.            phi_sav[chn][new][j] = phi[j];
  271. temp1=r[chn][new][j] * cos((double) phi[j]) - r_prime * cos((double) phi_prime);
  272. temp2=r[chn][new][j] * sin((double) phi[j]) - r_prime * sin((double) phi_prime);
  273.            temp3=r[chn][new][j] + fabs((double)r_prime);
  274.            if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;
  275.            else c[j] = 0;
  276.         }
  277. /*****************************************************************************
  278.  * Calculate the grouped, energy-weighted, unpredictability measure,         *
  279.  * grouped_c[], and the grouped energy. grouped_e[]                          *
  280.  *****************************************************************************/
  281.         for(j=1;j<CBANDS;j++){
  282.            grouped_e[j] = 0;
  283.            grouped_c[j] = 0;
  284.         }
  285.         grouped_e[0] = energy[0];
  286.         grouped_c[0] = energy[0]*c[0];
  287.         for(j=1;j<HBLKSIZE;j++){
  288.            grouped_e[partition[j]] += energy[j];
  289.            grouped_c[partition[j]] += energy[j]*c[j];
  290.         }
  291. /*****************************************************************************
  292.  * convolve the grouped energy-weighted unpredictability measure             *
  293.  * and the grouped energy with the spreading function, s[j][k]               *
  294.  *****************************************************************************/
  295.         for(j=0;j<CBANDS;j++){
  296.            ecb[j] = 0;
  297.            cb[j] = 0;
  298.            for(k=0;k<CBANDS;k++){
  299.               if(s[j][k] != 0.0){
  300.                  ecb[j] += s[j][k]*grouped_e[k];
  301.                  cb[j] += s[j][k]*grouped_c[k];
  302.               }
  303.            }
  304.            if(ecb[j] !=0)cb[j] = cb[j]/ecb[j];
  305.            else cb[j] = 0;
  306.         }
  307. /*****************************************************************************
  308.  * Calculate the required SNR for each of the frequency partitions           *
  309.  *         this whole section can be accomplished by a table lookup          *
  310.  *****************************************************************************/
  311.         for(j=0;j<CBANDS;j++){
  312.            if(cb[j]<.05)cb[j]=0.05;
  313.            else if(cb[j]>.5)cb[j]=0.5;
  314.            tb = -0.434294482*log((double) cb[j])-0.301029996;
  315.            bc[j] = tmn[j]*tb + nmt*(1.0-tb);
  316.            k = cbval[j] + 0.5;
  317.            bc[j] = (bc[j] > bmax[k]) ? bc[j] : bmax[k];
  318.            bc[j] = exp((double) -bc[j]*LN_TO_LOG10);
  319.         }
  320. /*****************************************************************************
  321.  * Calculate the permissible noise energy level in each of the frequency     *
  322.  * partitions. Include absolute threshold and pre-echo controls              *
  323.  *         this whole section can be accomplished by a table lookup          *
  324.  *****************************************************************************/
  325.         for(j=0;j<CBANDS;j++)
  326.            if(rnorm[j] && numlines[j])
  327.               nb[j] = ecb[j]*bc[j]/(rnorm[j]*numlines[j]);
  328.            else nb[j] = 0;
  329.         for(j=0;j<HBLKSIZE;j++){
  330. /*temp1 is the preliminary threshold */
  331.            temp1=nb[partition[j]];
  332.            temp1=(temp1>absthr[j])?temp1:absthr[j];
  333. /*do not use pre-echo control for layer 2 because it may do bad things to the*/
  334. /*  MUSICAM bit allocation algorithm                                         */
  335.            if(lay==1){
  336.               fthr[j] = (temp1 < lthr[chn][j]) ? temp1 : lthr[chn][j];
  337.               temp2 = temp1 * 0.00316;
  338.               fthr[j] = (temp2 > fthr[j]) ? temp2 : fthr[j];
  339.            }
  340.            else fthr[j] = temp1;
  341.            lthr[chn][j] = LXMIN*temp1;
  342.         }
  343. /*****************************************************************************
  344.  * Translate the 512 threshold values to the 32 filter bands of the coder    *
  345.  *****************************************************************************/
  346.         for(j=0;j<193;j += 16){
  347.            minthres = 60802371420160.0;
  348.            sum_energy = 0.0;
  349.            for(k=0;k<17;k++){
  350.               if(minthres>fthr[j+k])minthres = fthr[j+k];
  351.               sum_energy += energy[j+k];
  352.            }
  353.            snrtmp[i][j/16] = sum_energy/(minthres * 17.0);
  354.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  355.         }
  356.         for(j=208;j<(HBLKSIZE-1);j += 16){
  357.            minthres = 0.0;
  358.            sum_energy = 0.0;
  359.            for(k=0;k<17;k++){
  360.               minthres += fthr[j+k];
  361.               sum_energy += energy[j+k];
  362.            }
  363.            snrtmp[i][j/16] = sum_energy/minthres;
  364.            snrtmp[i][j/16] = 4.342944819 * log((double)snrtmp[i][j/16]);
  365.         }
  366. /*****************************************************************************
  367.  * End of Psychoacuostic calculation loop                                    *
  368.  *****************************************************************************/
  369.      }
  370.      for(i=0; i<32; i++){
  371.         if(lay==2)
  372.            snr32[i]=(snrtmp[0][i]>snrtmp[1][i])?snrtmp[0][i]:snrtmp[1][i];
  373.         else snr32[i]=snrtmp[0][i];
  374.      }
  375.      break;
  376.   case 3:
  377.      printf("layer 3 is not currently supported\n");
  378.      break;
  379.   default:
  380.      printf("error, invalid MPEG/audio coding layer: %d\n",lay);
  381.  }
  382.  
  383. /* These mem_free() calls must correspond with the mem_alloc() calls     */
  384. /* used at the beginning of this function to simulate "automatic"        */
  385. /* variables placed on the stack.                                        */
  386.  
  387.  mem_free((void **) &grouped_c);
  388.  mem_free((void **) &grouped_e);
  389.  mem_free((void **) &nb);
  390.  mem_free((void **) &cb);
  391.  mem_free((void **) &ecb);
  392.  mem_free((void **) &bc);
  393.  mem_free((void **) &wsamp_r);
  394.  mem_free((void **) &wsamp_i);
  395.  mem_free((void **) &phi);
  396.  mem_free((void **) &energy);
  397.  mem_free((void **) &c);
  398.  mem_free((void **) &fthr);
  399.  mem_free((void **) &snrtmp);
  400. }
  401.  
  402. /******************************************************************************
  403. routine to read in absthr table from a file.
  404. ******************************************************************************/
  405.  
  406. void read_absthr(absthr, table)
  407. FLOAT *absthr;
  408. int table;
  409. {
  410.  FILE *fp;
  411.  long j,index;
  412.  float a;
  413.  char t[80];
  414.  char ta[16];
  415.  
  416.  strcpy( ta, "absthr_0" );
  417.  
  418.  switch(table){
  419.     case 0 : ta[7] = '0';
  420.              break;
  421.     case 1 : ta[7] = '1';
  422.              break;
  423.     case 2 : ta[7] = '2';
  424.              break;
  425.     default : printf("absthr table: Not valid table number\n");
  426.  }
  427.  if(!(fp = OpenTableFile(ta) ) ){
  428.     printf("Please check %s table\n", ta);
  429.     exit(1);
  430.  }
  431.  fgets(t, 150, fp);
  432.  sscanf(t, "table %ld", &index);
  433.  if(index != table){
  434.     printf("error in absthr table %s",ta);
  435.     exit(1);
  436.  }
  437.  for(j=0; j<HBLKSIZE; j++){
  438.     fgets(t,80,fp);
  439.     sscanf(t,"%f", &a);
  440.     absthr[j] =  a;
  441.  }
  442.  fclose(fp);
  443. }
  444.